https://github.com/cran/VineCopula - Copula family
setwd('~/Documents/VaR')
prices <- read.csv('./data/StockIndexData.csv')
library(VineCopula)
library(copula)
##
## Attaching package: 'copula'
## The following object is masked from 'package:VineCopula':
##
## fitCopula
library(tseries)
prices.DJIA <- prices[,c('DJIA')]
summary(prices.DJIA)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6547 10220 11020 12020 13270 18640
n <- length(prices.DJIA)
rtn.DJIA <- rep(0,n-1)
for(i in seq(n-1)){
rtn.DJIA[i] <- log(prices.DJIA[i+1]/prices.DJIA[i]) * 100.
}
prices.GSPC <- prices[,c('GSPC')]
summary(prices.GSPC)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 676.5 1133.0 1296.0 1364.0 1478.0 2190.0
n <- length(prices.GSPC)
rtn.GSPC <- rep(0,n-1)
for(i in seq(n-1)){
rtn.GSPC[i] <- log(prices.GSPC[i+1]/prices.GSPC[i]) * 100.
}
prices.NDX <- prices[,c('NDX')]
summary(prices.NDX)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 938.5 1569.0 1944.0 2394.0 2954.0 4910.0 1063
n <- length(prices.NDX)
rtn.NDX <- rep(0,n-1)
for(i in seq(n-1)){
rtn.NDX[i] <- log(prices.NDX[i+1]/prices.NDX[i]) * 100.
}
prices.GDAXI <- prices[,c('GDAXI')]
summary(prices.GDAXI)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2203 4934 6172 6461 7583 12370 60
n <- length(prices.GDAXI)
rtn.GDAXI <- rep(0,n-1)
for(i in seq(n-1)){
rtn.GDAXI[i] <- log(prices.GDAXI[i+1]/prices.GDAXI[i]) * 100.
}
prices.FCHI <- prices[,c('FCHI')]
summary(prices.FCHI)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2403 3649 4220 4291 4836 6857 50
n <- length(prices.FCHI)
rtn.FCHI <- rep(0,n-1)
for(i in seq(n-1)){
rtn.FCHI[i] <- log(prices.FCHI[i+1]/prices.FCHI[i]) * 100.
}
prices.SSEC <- prices[,c('SSEC')]
summary(prices.SSEC)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1011 1572 2104 2288 2850 6092 317
n <- length(prices.SSEC)
rtn.SSEC <- rep(0,n-1)
for(i in seq(n-1)){
rtn.SSEC[i] <- log(prices.SSEC[i+1]/prices.SSEC[i]) * 100.
}
prices.SENSEX <- prices[,c('SENSEX')]
summary(prices.SENSEX)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2600 4757 13670 13040 18840 29680 202
n <- length(prices.SENSEX)
rtn.SENSEX <- rep(0,n-1)
for(i in seq(n-1)){
rtn.SENSEX[i] <- log(prices.SENSEX[i+1]/prices.SENSEX[i]) * 100.
}
rtn.DJIA.ar <- ar(rtn.DJIA, order.max = 1, method = 'ols')
rtn.DJIA.ar.resid.garch <- garch(rtn.DJIA.ar$resid[2:length(rtn.DJIA.ar$resid)],trace=FALSE)
rtn.GSPC.ar <- ar(rtn.GSPC, order.max = 1, method = 'ols')
rtn.GSPC.ar.resid.garch <- garch(rtn.GSPC.ar$resid[2:length(rtn.GSPC.ar$resid)],trace=FALSE)
rtn.NDX.1<-na.omit(rtn.NDX)
rtn.NDX.ar <- ar(rtn.NDX.1, order.max = 1, method = 'ols')
rtn.NDX.ar.resid.garch <- garch(rtn.NDX.ar$resid[2:length(rtn.NDX.ar$resid)],trace=FALSE)
rtn.GDAXI.1<-na.omit(rtn.GDAXI)
rtn.GDAXI.ar <- ar(rtn.GDAXI.1, order.max = 1, method = 'ols')
rtn.GDAXI.ar.resid.garch <- garch(rtn.GDAXI.ar$resid[2:length(rtn.GDAXI.ar$resid)],trace=FALSE)
rtn.FCHI.1<-na.omit(rtn.FCHI)
rtn.FCHI.ar <- ar(rtn.FCHI.1, order.max = 1, method = 'ols')
rtn.FCHI.ar.resid.garch <- garch(rtn.FCHI.ar$resid[2:length(rtn.FCHI.ar$resid)],trace=FALSE)
rtn.SSEC.1<-na.omit(rtn.SSEC)
rtn.SSEC.ar <- ar(rtn.SSEC.1, order.max = 1, method = 'ols')
rtn.SSEC.ar.resid.garch <- garch(rtn.SSEC.ar$resid[2:length(rtn.SSEC.ar$resid)],trace=FALSE)
rtn.SENSEX.1<-na.omit(rtn.SENSEX)
rtn.SENSEX.ar <- ar(rtn.SENSEX.1, order.max = 1, method = 'ols')
rtn.SENSEX.ar.resid.garch <- garch(rtn.SENSEX.ar$resid[2:length(rtn.SENSEX.ar$resid)],trace=FALSE)
rm(prices.DJIA, prices.GSPC, prices.NDX, prices.GDAXI, prices.FCHI,prices.SSEC, prices.SENSEX)
rm(rtn.DJIA, rtn.GSPC, rtn.NDX, rtn.GDAXI, rtn.FCHI,rtn.SSEC, rtn.SENSEX)
rm(rtn.NDX.1, rtn.GDAXI.1, rtn.FCHI.1,rtn.SSEC.1, rtn.SENSEX.1)
std.residual.DJIA <- rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.DJIA.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.DJIA.ar$resid/sqrt(rtn.DJIA.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GSPC <- rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GSPC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GSPC.ar$resid/sqrt(rtn.GSPC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.NDX <- rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.NDX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.NDX.ar$resid/sqrt(rtn.NDX.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.GDAXI <- rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.GDAXI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.GDAXI.ar$resid/sqrt(rtn.GDAXI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.FCHI <- rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.FCHI.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.FCHI.ar$resid/sqrt(rtn.FCHI.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SSEC <- rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SSEC.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SSEC.ar$resid/sqrt(rtn.SSEC.ar.resid.garch$fitted.values):
## longer object length is not a multiple of shorter object length
std.residual.SENSEX <- rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch$fitted.values)
## Warning in sqrt(rtn.SENSEX.ar.resid.garch$fitted.values): NaNs produced
## Warning in rtn.SENSEX.ar$resid/sqrt(rtn.SENSEX.ar.resid.garch
## $fitted.values): longer object length is not a multiple of shorter object
## length
rtns <- cbind(std.residual.DJIA[1067:4532],
std.residual.GSPC[1067:4532],
std.residual.NDX[2:3467],
std.residual.GDAXI[953:4418],
std.residual.FCHI[967:4432],
std.residual.SSEC[658:4123],
std.residual.SENSEX[674:4139]
)
cor(rtns,method='pearson')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00000000 0.97482107 0.4626361048 -0.02897670 -0.0448450319
## [2,] 0.97482107 1.00000000 0.4877918823 -0.02413000 -0.0454992144
## [3,] 0.46263610 0.48779188 1.0000000000 -0.05170464 0.0005284583
## [4,] -0.02897670 -0.02413000 -0.0517046389 1.00000000 -0.0142933609
## [5,] -0.04484503 -0.04549921 0.0005284583 -0.01429336 1.0000000000
## [6,] -0.02312062 -0.02512210 0.0025580560 0.01929367 0.0158593348
## [7,] -0.01925440 -0.02025158 -0.0048500626 0.01343810 -0.0123316220
## [,6] [,7]
## [1,] -0.023120619 -0.019254398
## [2,] -0.025122102 -0.020251582
## [3,] 0.002558056 -0.004850063
## [4,] 0.019293672 0.013438099
## [5,] 0.015859335 -0.012331622
## [6,] 1.000000000 0.013798025
## [7,] 0.013798025 1.000000000
cor(rtns,method='kendall')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.840936111 0.339267541 -0.018745030 -0.020497615
## [2,] 0.840936111 1.000000000 0.370060343 -0.012876103 -0.020819688
## [3,] 0.339267541 0.370060343 1.000000000 -0.026871801 -0.001039327
## [4,] -0.018745030 -0.012876103 -0.026871801 1.000000000 -0.012353858
## [5,] -0.020497615 -0.020819688 -0.001039327 -0.012353858 1.000000000
## [6,] -0.003618578 -0.004117842 -0.000861471 0.010674214 -0.002154094
## [7,] -0.004118841 -0.003936654 -0.004190783 -0.001609867 -0.020418012
## [,6] [,7]
## [1,] -0.003618578 -0.004118841
## [2,] -0.004117842 -0.003936654
## [3,] -0.000861471 -0.004190783
## [4,] 0.010674214 -0.001609867
## [5,] -0.002154094 -0.020418012
## [6,] 1.000000000 0.009381258
## [7,] 0.009381258 1.000000000
cor(rtns,method='spearman')
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.962107436 0.454675399 -0.027960277 -0.030258828
## [2,] 0.962107436 1.000000000 0.484046514 -0.019486288 -0.030636594
## [3,] 0.454675399 0.484046514 1.000000000 -0.040064311 -0.001568855
## [4,] -0.027960277 -0.019486288 -0.040064311 1.000000000 -0.018482587
## [5,] -0.030258828 -0.030636594 -0.001568855 -0.018482587 1.000000000
## [6,] -0.005895352 -0.006491070 -0.001236738 0.015962569 -0.002794678
## [7,] -0.006206004 -0.006164209 -0.006189422 -0.002249416 -0.030567593
## [,6] [,7]
## [1,] -0.005895352 -0.006206004
## [2,] -0.006491070 -0.006164209
## [3,] -0.001236738 -0.006189422
## [4,] 0.015962569 -0.002249416
## [5,] -0.002794678 -0.030567593
## [6,] 1.000000000 0.013779081
## [7,] 0.013779081 1.000000000
#pairs(rtns)
splom2(rtns)
g <- rtns[,1]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 1 10 70 366 1183 1401 377 53 3 2
##
## $density
## [1] 0.0002885170 0.0028851702 0.0201961916 0.1055972302 0.3413156376
## [6] 0.4042123485 0.1087709175 0.0152914022 0.0008655511 0.0005770340
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,2]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 4 14 82 362 1139 1409 398 52 4 2
##
## $density
## [1] 0.001154068 0.004039238 0.023658396 0.104443162 0.328620889
## [6] 0.406520485 0.114829775 0.015002885 0.001154068 0.000577034
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,3]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
##
## $counts
## [1] 1 3 13 132 423 1072 1259 464 90 7 1 1
##
## $density
## [1] 0.0002885170 0.0008655511 0.0037507213 0.0380842470 0.1220427005
## [6] 0.3092902481 0.3632429313 0.1338718984 0.0259665320 0.0020196192
## [11] 0.0002885170 0.0002885170
##
## $mids
## [1] -5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,4]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 3 19 119 431 1051 1258 478 88 16 3
##
## $density
## [1] 0.0008655511 0.0054818234 0.0343335257 0.1243508367 0.3032313907
## [6] 0.3629544143 0.1379111368 0.0253894980 0.0046162724 0.0008655511
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,5]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 2 20 122 422 1118 1205 479 82 10 6
##
## $density
## [1] 0.000577034 0.005770340 0.035199077 0.121754183 0.322562031
## [6] 0.347663012 0.138199654 0.023658396 0.002885170 0.001731102
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,6]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6
##
## $counts
## [1] 1 1 10 42 112 424 1099 1166 460 123 22 4 2
##
## $density
## [1] 0.000288517 0.000288517 0.002885170 0.012117715 0.032313907
## [6] 0.122331218 0.317080208 0.336410848 0.132717830 0.035487594
## [11] 0.006347374 0.001154068 0.000577034
##
## $mids
## [1] -6.5 -5.5 -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 5.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
g <- rtns[,7]
h<-hist(g, breaks=10, density=10, col="lightgray", xlab="Accuracy", main="Overall")
xfit<-seq(min(g),max(g),length=40)
yfit<-dnorm(xfit,mean=mean(g),sd=sd(g))
yfit <- yfit*diff(h$mids[1:2])*length(g)
lines(xfit, yfit, col="black", lwd=2)
h
## $breaks
## [1] -5 -4 -3 -2 -1 0 1 2 3 4 5
##
## $counts
## [1] 3 3 65 373 1228 1346 396 48 3 1
##
## $density
## [1] 0.0008655511 0.0008655511 0.0187536065 0.1076168494 0.3542989036
## [6] 0.3883439123 0.1142527409 0.0138488171 0.0008655511 0.0002885170
##
## $mids
## [1] -4.5 -3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 4.5
##
## $xname
## [1] "g"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
a <- rtns[,1]
b <- rtns[,2]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 7
## Name: BB1
##
## Parameter(s)
## ------------
## par: 1.05 (SE = 0.06)
## par2: 4.26 (SE = 0.1)
##
## Dependence measures
## -------------------
## Kendall's tau: 0.85 (empirical = 0.84, p value < 0.01)
## Upper TD: 0.82
## Lower TD: 0.86
##
## Fit statistics
## --------------
## logLik: 5119.14
## AIC: -10234.27
## BIC: -10221.97
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## 0.9714474
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## 0.9705497 4.1083364
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df = df),
list(df = df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## 7.36012
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2),margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
a <- rtns[,1]
b <- rtns[,3]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 2
## Name: t
##
## Parameter(s)
## ------------
## par: 0.54 (SE = 0.01)
## par2: 2 (SE = 0.11)
##
## Dependence measures
## -------------------
## Kendall's tau: 0.37 (empirical = 0.34, p value < 0.01)
## Upper TD: 0.42
## Lower TD: 0.42
##
## Fit statistics
## --------------
## logLik: 808.31
## AIC: -1612.62
## BIC: -1600.31
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## 0.4660699
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## 0.5363481 1.7946408
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df=df),
list(df=df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## 0.836929
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
a <- rtns[,1]
b <- rtns[,4]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 2
## Name: t
##
## Parameter(s)
## ------------
## par: -0.03 (SE = 0.02)
## par2: 13.28 (SE = 3.49)
##
## Dependence measures
## -------------------
## Kendall's tau: -0.02 (empirical = -0.02, p value = 0.1)
## Upper TD: 0
## Lower TD: 0
##
## Fit statistics
## --------------
## logLik: 9.75
## AIC: -15.5
## BIC: -3.19
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## -0.03054467
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## -0.0300592 13.2810415
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df=df),
list(df=df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## -0.01112868
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
a <- rtns[,1]
b <- rtns[,5]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 2
## Name: t
##
## Parameter(s)
## ------------
## par: -0.03 (SE = 0.02)
## par2: 9.21 (SE = 1.66)
##
## Dependence measures
## -------------------
## Kendall's tau: -0.02 (empirical = -0.02, p value = 0.07)
## Upper TD: 0.01
## Lower TD: 0.01
##
## Fit statistics
## --------------
## logLik: 21.35
## AIC: -38.7
## BIC: -26.4
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## -0.03932221
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## -0.03382907 9.20959786
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df=df),
list(df=df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## -0.01643638
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
a <- rtns[,1]
b <- rtns[,6]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 2
## Name: t
##
## Parameter(s)
## ------------
## par: -0.01 (SE = 0.02)
## par2: 16.17 (SE = 5.01)
##
## Dependence measures
## -------------------
## Kendall's tau: -0.01 (empirical = 0, p value = 0.75)
## Upper TD: 0
## Lower TD: 0
##
## Fit statistics
## --------------
## logLik: 6.34
## AIC: -8.67
## BIC: 3.63
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## -0.01696391
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## -0.0129055 16.1649914
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df=df),
list(df=df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## -0.01336535
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
a <- rtns[,1]
b <- rtns[,7]
m <- pobs(as.matrix(cbind(a,b)))
u <- m[,1]
v <- m[,2]
mu.1 <- mean(a)
sd.1 <- sd(a)
mu.2 <- mean(b)
sd.2 <- sd(b)
plot(a,b,pch='.')
abline(lm(a~b),col='red',lwd=1)
selectedCopula <- BiCopSelect(u,v,familyset=NA)
summary(selectedCopula)
## Family
## ------
## No: 124
## Name: Rotated Tawn type 1 90 degrees
##
## Parameter(s)
## ------------
## par: -1.1 (SE = 0.08)
## par2: 0.06 (SE = 0.05)
##
## Dependence measures
## -------------------
## Kendall's tau: -0.01 (empirical = 0, p value = 0.72)
## Upper TD: 0
## Lower TD: 0
##
## Fit statistics
## --------------
## logLik: 4.14
## AIC: -4.28
## BIC: 8.02
n.cop <- normalCopula(dim=2)
set.seed(500)
fit <- fitCopula(n.cop,m,method='ml')
coef(fit)
## rho.1
## -0.01430268
rho <- coef(fit)[1]
copula_dist <- mvdc(copula=normalCopula(rho,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
t.cop <- tCopula(dim=2)
set.seed(500)
fit <- fitCopula(t.cop,m,method='ml')
coef(fit)
## rho.1 df
## -0.01250382 34.12112734
rho <- coef(fit)[1]
df <- coef(fit)[2]
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("t","t"),
paramMargins=list(list(df=df),
list(df=df)))
sim1 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
c.cop <- claytonCopula(dim=2)
set.seed(500)
fit <- fitCopula(c.cop,m,method='ml')
coef(fit)
## param
## -0.01099794
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=claytonCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim2 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
g.cop <- gumbelCopula(dim=2)
set.seed(500)
m <- pobs(as.matrix(cbind(rtns[,1],rtns[,3])))
fit <- fitCopula(g.cop,m,method='ml')
coef(fit)
## param
## 1.521606
theta <- coef(fit)[1]
copula_dist <- mvdc(copula=gumbelCopula(theta,dim=2), margins=c("norm","norm"),
paramMargins=list(list(mean=mu.1, sd=sd.1),
list(mean=mu.2, sd=sd.2)))
sim3 <- rmvdc(copula_dist, 3466)
## Warning: 'rmvdc' is deprecated.
## Use 'rMvdc' instead.
## See help("Deprecated")
plot(a,b,main='normal-copula')
points(sim[,1],sim[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='t-copula')
points(sim1[,1],sim1[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Clayton copula')
points(sim2[,1],sim2[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)
plot(a,b,main='Gumbel copula')
points(sim3[,1],sim3[,2],col='red')
legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)